Mon. Not. R. Astron. Soc. 000, Q-?? (0000) Printed 1 February 2008 (MN WT^X style file vl.4) 



Hydrodynamic Simulations of Propagating Warps and 
Bending Waves in Accretion Discs 

Richard P. Nelson * & John C.B. Papaloizou t 

Astronomy Unit, Queen Mary & Westfield College, Mile End Rd, London El J^NS 



Received *****; in original form 



ABSTRACT 

We present the results of a study of propagating warp or bending waves in accretion 
discs. Three dimensional hydrodynamic simulations were performed using SPH, and 
the results are compared with calculations based on the linear theory of warped discs. 
We examine the response of a gaseous disc to an initially imposed warping disturbance 
under a variety of physical conditions. We consider primarily the physical regime in 
which the dimensionless viscosity parameter a < H/r, where H/r is the disc aspect 
ratio, so that bending waves are expected to propagate. We also performed calculations 
for disc models in which a > H/r, where the warps are expected to evolve diffusively. 
Small amplitude (linear) perturbations are studied in both Keplerian and slightly non 
Keplerian discs, and we find that the results of the SPH calculations can be reasonably 
well fitted by those of the linear theory. The main results of these calculations are: 
(i) the warp in Keplerian discs when a < H/r propagates with little dispersion, 
and damps at a rate expected from estimates of the code viscosity, {ii} warps evolve 
diffusively when a > H/r, {in) the slightly non Keplerian discs lead to a substantially 
more dispersive behaviour of the warps, which damp at a similar rate to the Keplerian 
case, when a < H/r. 

Initially imposed higher amplitude, nonlinear warping disturbances were studied in 
Keplerian discs. The results indicate that nonlinear warps can lead to the formation 
of shocks, and that the evolution of the warp becomes less wave-like and more diffusive 
in character. 

This work is relevant to the study of the warped accretion discs that may occur around 
Kerr black holes or in misaligned binary systems, and is mainly concerned with discs 
in which a < H/r. The results indicate that SPH can model the hydrodynamics of 
warped discs, even when using rather modest numbers of particles. 

Key words: accretion discs-bending waves-warps-SPH-simulations 



1 INTRODUCTION 

Accretion discs occur in a variety of astrophysical con- 
texts such as around T-Tauri stars, around Galactic X-ray 
sources such as Cygnus X-1, or around massive black holes 
at the centres of active galactic nuclei. In a situation where 
these discs are subject to out-of-plane forcing, the disc may 
become globally twisted or warped. Examples of where this 
may occur are in a binary system where the binary or- 
bit plane is misaligned with respect to the disc midplane 
(e.g. Papaloizou & Terquem 1995, Larwood et al. 1996), or 
around a rotating black hole where the Lense-Thirring ef- 



* Email address: R.P.Nelson@qmw.ac.uk 

t Email address; J.C.B.Papaloizou@qmw.ac.uk 



feet may cause a warped disc to form (Bardeen & Petterson 
1975, Kumar & Pringle 1985, Nelson & Papaloizou 1999). 

The linear theory of warped discs, in which both vis- 
cous and pressure effects were included, was studied by Pa- 
paloizou & Pringle (1983). They examined the evolution of 
discs in the regime where a > H/r, a being the dimen- 
sionless viscosity coefficient of Shakura & Sunyaev (1973) 
and H/r being the disc aspect ratio, and showed that warps 
evolve diffusively in this regime (see also Ogilvie 1999). The 
diffusion coefficient associated with the warp evolution was 
found to be larger than that associated with mass flow 
through the disc by a factor of ~ l/(2a^), assuming an 
isotropic viscosity, indicating that disc warps diffuse more 
rapidly than had previously been supposed (e.g. Bardeen 
& Petterson 1975). For discs in which a < H/r, the gov- 
erning equation for disc warp evolution changes from be- 
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ing a diffusion equation to a wave equation, indicating that 
warps propagate as bending waves in this physical regime. 
These bending waves were studied using hnear theory by 
Papaloizou & Lin (1994, 1995), who showed under the as- 
sumption of very low frequency waves that warping distur- 
bances in a Keplerian disc propagate with no dispersion at a 
speed corresponding to an appropriate average of the sound 
speed. 

The dynamics of warped discs have recently been stud- 
ied using numerical simulations performed with SPH. Lar- 
wood et al. (1996) examined the structure of a disc around 
one component of a binary system in which the disc mid- 
plane and the binary orbit plane were misaligned. A similar 
study of warped circumbinary discs was undertaken by Lar- 
wood & Papaloizou (1997). In a companion paper to this 
one (Nelson & Papaloizou 1999), we study the dynamics of 
warped accretion discs orbiting around rotating black holes, 
where the angular momentum vector of the black hole is 
misaligned with that of the disc. In this case, the combined 
effects of differential Lense-Thirring precession due to the 
Kerr black hole and viscosity cause the inner regions of these 
discs to become warped. 

In this paper, we examine the evolution of warping dis- 
turbances in accretion discs using SPH in greater detail, and 
compare the results of the nonlinear simulations with solu- 
tions to the linearized initial value problem obtained using 
a standard finite difference method. The wavelengths of the 
warping perturbations that we study are on the order of 
the disc thickness, so we relax the assumption of low fre- 
quency bending waves in the linear theory. We examine the 
evolution of warps under a variety of physical conditions, 
including discs orbiting in both Keplerian and slightly non 
Keplerian central potentials, and warps with both small and 
large amplitudes. We find good agreement between the SPH 
simulations and the predictions of linear theory for low am- 
plitude warps. 

In the regime where a < H/r, the warps propagate 
across the disc as inward and outward moving bending 
waves. When a > H/r, the warp evolution becomes diffu- 
sive, as expected. We find that the evolution of large ampli- 
tude warps becomes increasingly diffusive and less wave-like 
due to nonlinear effects such as shocks, which form because 
perturbed velocities on the order of the sound speed are 
generated. 

The plan of this paper is as follows. In section (Q) we 
present the basic equations of the problem. In section (0) we 
present a discussion of the linear theory of gaseous warped 
discs. In section (^ we describe the numerical method used 
to perform the nonlinear simulations , and in section (^ 
we present the results of these simulations and those ob- 
tained from appropriate solutions of the linearized initial 
value problem. Finally, a discussion of our results and our 
conclusions are presented in section (tel). 



2 EQUATIONS OF MOTION 

In order to describe a compressible fluid, we adopt the con- 
tinuity and momentum equations in the form 



-VP -V'i> + S„,, 



1 + ^^-^ 
and 



(2) 



(3) 



H ' 
where 

d d ^ 
— = h v.V 

dt dt 

denotes the convective derivative, p is the density, v is the 
velocity, P is the pressure, "1> is the gravitational potential 
and S„isc represents the viscous force per unit mass. 

To describe the disc, we adopt a cylindrical coordinate 
system (r, 0, z) based on the central mass around which the 
accretion disc orbits. The corresponding Cartesian coordi- 
nates are {x, y, z). The distance to the central mass is given 
hy R = Vr^ + z^, and the position vector measured from 
there is denoted by r. 

The gravitational potential <E> used in the calculations 
is given by 



$ : 



GM 



(4) 



where h is an optional softernng parameter used to prevent 
numerical divergences in the gravitational force for accre- 
tion disc models that extend to the centre of the coordinate 
system. 

The equation of state is taken to be that of a polytrope 
with 7 = 5/3: 



(5) 







(1) 



Energy dissipated through the action of artificial viscosity 
is simply allowed to leave the system, so that a barotropic 
equation of state is assumed throughout. 



3 LINEAR THEORY OF WARPS AND 
BENDING WAVES 

The linear theory of warps in a vertically stratified gaseous 
accretion disc, including hydrodynamical effects, was ini- 
tially investigated by Papaloizou & Pringle (1983). 

The analysis presented in their paper, however, applied 
only to viscous accretion discs in which the dimensionless 
Shakura & Sunyaev (1973) viscosity parameter, a, was con- 
strained to lie in the range Hjr < a < 1, where H/r is the 
ratio of disc height to radius. Under these circumstances, 
globally warped discs were found to evolve diffusively on a 
time scale tr, ~ ci'tv, where — /v is the viscous time 
scale appropriate to accretion, expected from standard thin 
disc theory. We comment here that in the above a standard 
isotropic viscous stress tensor was assumed. 

For the regime a < H/r, however, the nature of the 
governing equation for the disc tilt changes from being of 
the diffusion type to being a wave equation, such that warps 
are communicated through the disc via the propagation of 
bending waves (Papaloizou & Lin 1994, 1995). 



3.1 Initial Value Problem 

Small amplitude warps may be considered to be linear per- 
turbations of a disc with midplane initially coincident with 
the {x, y) plane. The dependence of all perturbations may 
be taken into account through a factor e''™''*', where m is 
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the azimuthal mode number. For the global warps consid- 
ered here, m = 1. The components of the linearized equa- 
tions of motion for a barotropic fluid (see eg. Papaloizou & 
Lin 1994, 1995) may be written, denoting perturbations by 
a prime: 



dv' 



dr 
-f iQ,v' 



iW 

r 
dW 
dz 



(6) 



dt 

P'/p ~ p'cg/p, with Cs being the sound speed 



Here W 

The linearized continuity equation may be written as 



V dt 



dr 



dz 



(7) 



For a disc in a spherically symmetric external potential, 
because there is no preferred direction for the disc rotation 
axis, equations (^) and (Q) have a solution corresponding to 
a time independent rigid tilt. Then W — —irzQ,^g, where 
g is a constant inclination. The corresponding components 
of the velocity perturbation are v'^ = rflg,v'^ = —zflg,v'^ — 
—izg{d{rQ) / dr) . Here we recall that for a barotropic disc fl 
is a function of r only. 

For perturbations corresponding to large scale warps, 
we expect the local inclination to vary on a radial length 
scale significantly greater than the thickness H. Then it is 
reasonable to assume that it is also approximately indepen- 
dent of z. We remark that this was found to be the case 
in the linear calculations of bending waves by Papaloizou & 
Lin (1995), which took the vertical structure of the disc fully 
into account, even when the radial wavelength was compa- 
rable to the vertical thickness. Thus to describe bending 
waves with long radial wavelength, we set W = —irzQ^g, 
and assume the local inclination g and thus are approx- 
imately independent of z. Then (H) implies that v'^, and v'^ 
are both proportional to z. Therefore we set v'^ — zq'^, and 
v'^ = zq'^ respectively, where q'r, and q'^ are also assumed 
approximately independent of z. Equation (Mj then gives 



dq'r ^ .„ , „„ , . drn^g 

+ iQ.q^ + q^r ' = -Q g 

+ iflv'^ = irQ,^g (8) 

The continuity equation (^) is then multiplied by z and 
vertically integrated through the disc with the result 



IrQ- 
where 



1= / ^dz, 

Cs 



i d(rpq'^) pq'^ , 

^ \ hiEu^, (9) 

r or r 



(10) 



and 



P ~ pz dz. 



(11) 



Equations (|8|) and (|9|) provide a set of linear equations 
for the inclination, g, as a function of r and t. They may 
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be used to follow the evolution of a disturbance from initial 
data. We compare the time dependent evolution calculated 
in this way with that obtained from non-linear SPH calcu- 
lations below. 



3.2 Low Frequency Bending Waves 

For the purpose of discussing low frequency waves, it is con- 
venient to consider solutions of the linearized equations for 
which the t dependence is separated through a factor e^"^, 
where a is the mode frequency. 

The cornponents of the velocity perturbations are then 
given from (H) by (e.g. Papaloizou & Terquem 1995): 



z 



-{a + n)g- 



{rQ^{a + n)^+ga\3n + a)) 



(ct + f7)2 - «:2 



_ igQ +i{v'r/z)r 

z 



-1 d(r^n) 
dr 



n^r 



dz 



where 



2n d(r^n) 



(12) 



(13) 



r dr 

is the square of the epicyclic frequency. 

For low frequency modes in a near Keplerian disc with 
\a\ <C fl, it should be noted that there is a near-resonance 
through the near vanishing of {a + Q)^ — , with the result 
that slowly varying warps may induce large horizontal mo- 
tions and vertical shear -gf-'^ ■ When considering the 

effects of viscous dissipation, provided the viscosity is not 
highly anisotropic, it is the damping of this vertical shear 
that provides the dominant effect of viscosity. 

The inclination, g, is governed by a single, second order 
ordinary differential equation (e.g. Papaloizou & Lin 1994, 
Papaloizou & Terquem 1995). For the low frequency limit 
we are interested in, this may be approximated as: 

d ( pQ, dg^ 



4gn{ujz + a)I+ — (^ 



a + Q, — K dr 



= 0. 



(14) 



Here the free particle nodal precession frequency, luz , is given 
by 



2nL0z = ( - ^ * 



dz^ 



(15) 



In a near Keplerian disc, luz and the apsidal precession 
frequency SI — k are small compared to Q and low frequency 
disturbances may be considered for which a is of comparable 
magnitude. 

For a strictly Keplerian disc where ui^ ~ and k = fl, 
from (|l|) it may be seen that localised disturbances, for 
which g oc e''^'^ obey the dispersion relation given by 



2 2 7 2 

a = a K 
where 

1 [P 

2 V J 



(16) 



(17) 



It may be seen from the definition of p and I in equations 
(0) and (0) that is a mean sound speed, such that in 
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the limit ct — > 0, warping disturbances propagate without 
dispersion at half of this mean sound speed. For finite a, 
dispersive corrections are of order a/Q. 

The full set of equations and (^), which allow for per- 
turbations of arbitrary frequency, are used to calculate the 
evolution of linear bending waves below. However, the type 
of disturbances considered can in fact be well described using 
the low frequency approximation (see Nelson, Papaloizou & 
Terquem 1999) and so may be regarded as effectively decou- 
pled from other high frequency modes. 



3.3 Effect of a Small Viscosity 

In a Keplerian disc, the most important effect of a small 
viscosity is to act on the resonantly induced horizontal mo- 
tions through the (r, z) and {<j>, z) components of the viscous 
tensor. If the kinematic viscosity is taken to be 

P 



(18) 



where F is an arbitrary function of radius and Q = pzdz, 
then the main effect of the viscosity is to replace the resonant 
denominator 



in equation ( p^ ) by 



In a thin Keplerian disc. 



pzdz = P/Q 



Setting r = aiQ, we find u — {ctiP) / (pfi) which gives a pre- 
scription equivalent to that of Shakura & Sunyaev (1973). 
This scheme for incorporating viscosity is included in the 
governing equations (g) and (^ for the time dependent evo- 
lution of warping disturbances by the straightforward oper- 
ator replacement 



in equations (^) while leaving (^) unaltered. 

However, it should be noted that the ai here differs 
from the 'a' parameter usually considered in thin disc theory 
since it acts through the (r, z) and {(j>, z) components of the 
viscous tensor rather than through the (r, (j>) component. 
For this reason we introduce the subscript 1 from now on. 

When viscosity is introduced in this way, equation ( |l4| ) 
is modified to read 



4:gn{ujz 



ia\Q. + Q, — K dr 



= 



(19) 



The associated local dispersion relation (for a;z = Q. — k, = Q) 
is given through 



(20) 



o" — iaiQ, 



If |a| << QiSl in the low frequency limit, then equation (|20| 
becomes 

• 2 I 2 

la k 



This is the form of the dispersion relation one obtains from 
a diffusion equation which implies that for a Keplerian disc, 
warps evolve diffusively with diffusion coefficient given by 



V = 



401^ 



(21) 



The corresponding time scale for global evolution of the disc 
is given by 



r 

V' 



(22) 



The transition between wave-like and diffusive regimes 
is expected to occur when the diffusion time above, calcu- 
lated for the whole disc, becomes comparable to the wave 
propagation time across it, or when qi ~ H/r (e.g. Pa- 
paloizou & Lin 1994). 

It is very important to note that provided qi is signif- 
icantly below unity, in either regime the communication of 
warp information across the disc can be very much more 
rapid than the time scale for matter to accrete through it. 
This is in fact confirmed in our simulations and expected 
from Papaloizou & Pringle (1983). 

The non linear numerical simulations of polytropic discs 
presented in this paper apply for the most part to the pa- 
rameter regime in which qi < H/r, so that the time depen- 
dent behaviour of small warping disturbances should exhibit 
wave-like behaviour. To check that this is indeed the case, 
we have performed a number of calculations of the evolution 
of small amplitude bending waves in section (^. Results ob- 
tained with SPH are compared directly with the predictions 
of linear theory given by integrating equations (^ and (^ 
forward in time, starting from the appropriate initial condi- 
tions, and modified to includ e th e effects of viscosity in the 
manner described in section (3.3). 



4 NUMERICAL METHOD 

The set of fiuid equations described in section (^) are solved 
using smoothed particle hydrodynamics (Lucy 1977; Gin- 
gold & Monaghan 1977). SPH uses particles to represent 
a subset of the fiuid elements that arise in the Lagrangian 
description of a fluid. The version of SPH used in the cal- 
culations presented here is a conservative formulation of the 
method that employs variable smoothing lengths (Nelson 
& Papaloizou 1994). For the sake of brevity, we do not pro- 
vide a detailed description of the method, but only highlight 
those points salient to the work presented in this paper. A 
detailed description of the code is presented in Nelson & 
Papaloizou (1994), along with a number of test calculations. 

In order for the method to be conservative, the smooth- 
ing lengths must be functions only of the interparticle sepa- 
rations. Accordingly, we flnd the Ntol nearest neighbours 
and calculate the smoothing lengths, hi, using the expres- 



Nfar 
ri — 1 



(23) 



where the summation is over the Nfar most distant nearest 
neighbours of particle i. For the calculations presented here, 
we take Ntol = 45 and Nfar = 6. 

The numerical method employs an artificial viscosity 
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Figure 1. This figure sliows the variation of disc tilt versus radius, as a function of time, for a linear calculation {dotted line), and a 
non linear SPH calculation [solid line). The assumed parameters for the linear calculation were M = 8.33, ai = 0.04/ri/2. The SPH 
calculation is described in the text. 



term which allows shocks to be properly modeled. The mag- 
nitude of this viscosity is controlled by two parameters, 
asPH and /3s ph as defined by Nelson & Papaloizou (1994). 
We set these parameters to have the values asPH = 0.5 and 
PSPH = 0.0. 

It should be remarked that in addition to providing a 
bulk viscosity, the artificial viscosity provides a shear vis- 
cosity arising from the finite resolution of the method which 
allows particles to communicate across shearing interfaces. 
The particles also develop random velocities that result from 
stochastic fluctuations in their pressure forces, and these mo- 
tions lead to a diffusive evolution of the particle distribution. 
Test calculations that used ring spreading simulations (Lar- 
wood et al 1996), and more recent calibration runs (Bryden 
et al. 1999) indicate that these processes lead to the viscous 



evolution of accretion disc models calculated with the SPH 
code. The code viscosity is characterised by a range of val- 
ues for the dimensionless shear viscosity parameter, a, that 
depend on the value of H/r in the disc. The dimensionless 
parameter a is defined through the expression for the kine- 
matic viscosity v — aH^Q, where H/r = , with M 
being the midplane Mach number in the disc. This does not 
vary greatly over the radial extent of most of the models 
considered here. For disc models with = 10 (i.e. with as- 
pect ratio H/r ~ 0.1) we find that a ~ 0.02, and for models 
with M = 30 we find that a ~ 0.03 - 0.04. 
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Figure 2. This figure shows the variation of disc tilt versus radius, as a function of time, for a linear calculation (dotted line), and a non 
linear SPH calculation {solid line). The assumed parameters for the linear calculation were M = 10, ai = 0.0. The SPH calculation is 
described in the text. 



5 BENDING WAVE CALCULATIONS 

We have performed a number of calculations in which small 
amplitude, localised warping disturbances were applied to 
accretion discs, and then allowed to evolve. The results of 
these calculations were compared with the results obtained 
from the linear theory through equations (^) and (^) dis- 
cussed in section These simulations, and their compar- 
ison with linear calculations which used a surface density 
profile matching that of the discs, are presented below. 

5.1 Initial and Boundary Conditions 

The calculations presented in this section employed differ- 
ing numbers of particles, and were performed for two dif- 
ferent values of the gravitational softening parameter, b, in- 



troduced in section The central mass in each of these 
calculations was taken to be M = 1, with the gravitational 
constant G — 1. The initial conditions for each calculation 
performed are described below. 

A number of calculations were performed for discs or- 
biting in purely Keplerian potentials, in which the soften- 
ing parameter 6 = 0. In these cases, the disc models were 
constructed by creating a gaseous annulus, with a free in- 
ner boundary at r = rin = 0.2 and a free outer boundary 
at r = Tout = 0.8, rotating about the z axis. This annu- 
lus contained A'^ — 102, 000 particles, and the required disc 
thickness (or equivalently disc midplane Mach number) was 
obtained by setting the polytropic constant, K, to the ap- 
propriate value. These disc models were then allowed to re- 
lax by evolving them for a time that was slightly greater 
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Figure 3. This figure sliows the variation of disc tilt versus radius, as a function of time, for a hnear calculation {dotted line), and a 
non linear SPH calculation {solid line). The assumed parameters for the linear calculation were M = 10, ai = 0.04/ri/2. The SPH 
calculation is described in the text. 



than two orbital periods at their outside edge. The surface 
density was initially constant, and was found to evolve very 
little during this process of relaxation. 

The bending wave calculations for these disc models 
were initiated by tilting an annulus in the disc, which lay 
between r = rmin = 0.4 and r = Tmax = 0.6. The tilt was 
imposed by rotating about the diameter coincident with the 
X axis. The tilt amplitude varied in radius according to the 
expression 



g(r) = go sin 



7r(r ■ 



Ar 



(24) 



where Ar = Tmax — fmin- Calculations were performed for 
discs in which the midplane Mach number was TVf = 10 and 
go took the values of 5, 10, and 20 degrees. An additional cal- 



culation was performed in which the midplane Mach number 
= 30 and go = ^ degrees. 

A suite of calculations were performed in which the discs 
extended all the way to i? = 0, so that softening of the 
central potential was required. In these cases, units were 
chosen such that the gravitational softening parameter b = 
0.2, and the disc outer radius R = 1. Calculations in this 
case used differing number of particles, with A'' — 20000, 
40000, and 102, 000. The disc thickness (or mid plane Mach 
number) was obtained by setting an appropriate value of K, 
and the initial disc tilt was obtained using equation (p^, 
with Vmin ~ 0.5 and rmax ~ 0.8. Tilt amplitudes with go = 
3, 7, 10, and 15 degrees were considered, with the midplane 
Mach number being A4 = 10. When discussing the results 
of these calculations below, we will concentrate primarily on 



© 0000 RAS, MNRAS 000, [|-?? 



8 Richard P. Nelson and John C.B. Papaloizou 




-0.7 -0.6 -0.5 -0.4 -0.3 -0.7 -0.6 -0.5 -0.4 -0.3 




Figure 4. This figure shows the perturbed radial velocity arising from bending wave propagation through the disc. Note the vertical 
shear as perturbed radial motion, which is odd in z, is generated in the disc by the warp. The arrows at the bottom of each panel indicate 
the magnitude of the midplane sound speed at each radius. 



the cases with go = 7 degrees, and N = 20000 and 102, 000 
particles, in order to compare the effects of changing the 
particle number and resolution on the results. 



primarily on a low amplitude tilt and a high amplitude tilt 
calculation. 



5.2 Results 

We now present and compare the results from a subset of the 
free bending wave calculations. We first discuss the results 
for a low amplitude (go = 5°) warping disturbance applied to 
a Keplerian disc. We also present results for two calculations 
performed with discs orbiting in softened potentials, but in 
which the particle numbers used differed from one another. 
Finally, we compare the results of three calculations in which 
the initial tilt amplitudes differed from one another, focusing 



5.2.1 Low Amplitude Tilt in Keplerian Disc 

The plots presented in Fig. (1) show the time evolution of 
the disc tilt g{r) during a calculation performed for a disc 
with M = 10 and N — 102000 particles. When discussing 
the figures in this paper, we will use the convention that 
the top left panel will be referred to as panel 1, with the 
remaining panels being labelled as 2, 3, and 4 when moving 
from left to right and from top to bottom. The time corre- 
sponding to each panel is shown in the top right hand corner 
in units of Q,~^ evaluated at i? = 1. The initial tilt applied 
to the disc was of the form given by equation (^41), with 
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Figure 5. This figure sliows tiie radial variation of tilt angle, as a function of time, for a bending wave calculation performed for a disc 
orbiting in a softened potential. Tlie solid line shows the results from the non linear SPH calculation, which is described in the text. The 
dot-dashed line shows the evolution of a linear calculation, with assumed parameters Ai = 10, ai = 



go — 5°, and may be seen in panel 1 of Fig. (1). The solid 
line in Fig. (1) represents the variation of tilt angle g with 
R for the non linear SPH calculation, whereas the dashed 
line represents a linear calculation computed using the time 
dependent linearized equations described in section 

Moving from left to right and from top to bottom, we 
can observe that the initial 'pulse' starts to broaden as the 
locally applied disc warp begins to propagate in towards the 
disc centre and out towards its exterior edge. After a time of 
t = 2.33, the initial pulse can be seen to be splitting into two 
separate pulses, corresponding to an in-going and an out- 
going bending wave, as one would expect from linear theory. 
After a time oit = 3.47 the waves have more or less reached 
the outer edge of the disc, and are almost fully separated. 
The parameters used in fitting the linear calculation to the 



non linear SPH calculation assumed a disc model in which 
the midplane Mach number A4 = 8.33, and took a vari- 
able value of the viscosity parameter ai = O.OA/r^^^, where 
this radial dependence of the viscosity arises from assuming 
an isotropic viscosity, and a steady state disc with uniform 
surface density E, aspect ratio (H/r), and mass flux. 

In addition to the fit presented in Fig. (1), we also pro- 
vide a comparison between the nonlinear SPH calculation 
and two further linear calculations in Figs. (2) and (3). The 
linear calculation in Fig. (2) assumed a disc midplane Mach 
number A4 = 10, and neglected the effects of viscosity en- 
tirely (i.e ai = 0.0). It is apparent that the SPH calculation 
mimics the broad features of the linear calculation in a gen- 
eral sense, but that a good fit is not obtained since viscosity 
is required to damp the amplitude of the waves. The lin- 
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Figure 6. This figure sliows tlie radial variation of tilt angle, as a function of time, for two bending wave calculations performed with 
discs orbiting in softened potentials. The solid line shows the results from a non linear SPH calculation using 20,000 particles, the dashed 
line shows results from a calculation employing 102,000 particles. 



ear calculation presented in Fig. (3) assumed = 10 and 
Qi = In this case it appears that the inward prop- 

agating pulse produced by the SPH calculation travels too 
fast, and the leading edge is not particularly well fitted by 
the linear calculation presented. The outward moving pulses 
in both the linear and nonlinear runs appear to match very 
closely, however, as shown in panel 4 of Fig. (3). We remark 
that we expect the results of the SPH calculations to be 
more accurate for the outward moving bending waves since 
the resolution is better in the outer parts of the disc models 
than in the inner parts. 

All of the linear calculations presented in Figs. (1) - 
(3) provide reasonable fits to the nonlinear SPH calcula- 
tion. The region in which the largest discrepancy arises in 
the Figs. (1) and (2) cases is at a radius of 71 ~ 0.5 - i.e. 



the radial position about which the initial warping distur- 
bance was centred. The linear calculations predict that the 
tilt angle, g, should drop sharply in this region as the in- 
going and out-going waves separate, with the magnitude of 
this decrease in g being determined by Oi. A larger viscosity 
leads to a more diffusive evolution of the warp, and thus to 
pulses that separate less cleanly. It is possible that the effec- 
tive viscosity (i.e. qi) in the SPH calculations is a function 
of the local conditions, and that the particle relaxation in 
the region initially occupied by the applied warping distur- 
bance leads to a larger effective viscosity there which is not 
accounted for in the simple steady state form of the viscosity 
adopted in the linear calculations. In our experience a very 
precise fit to the simulations may be obtained if an addi- 
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tional radial dependence is allowed in ai and A4. However, 
we do not pursue this here. 

The form of the perturbed velocity field expected to be 
induced in the disc when it is warped is given by equation 
(^). It may be seen that as a bending wave travels through 
the disc, it leads to the excitation of radial motions that 
are odd functions of z, such that a vertical shear is induced. 
This shearing motion may be observed in Fig. (4). Each of 
the panels shows particle projections onto the y-z plane for 
a thin cross section of the disc centred on the y axis. The 
disc is shown for —0.8 < y < —0.2. The radial velocities 
(vr) of the particles are indicated by an arrow, but only for 
particles for which \vr\ is above a small threshold value. The 
length of the arrow indicates the magnitude of the velocity, 
and the arrows located at the bottom of each panel show 
the magnitude of the midplane sound speed at each radial 
position in the disc. The top panel shows the initial state 
of the disc, with the applied warping disturbance centred 
on y — —0.5 being apparent. The residual velocities in this 
panel are just due to the residual random motions of the par- 
ticles. The middle panel shows the disc after it has evolved 
for t — 1.14, and it is obvious that vertical shearing motions 
have been excited by the initial warp. After further evolu- 
tion, the initial pulse separates into an inward and outward 
travelling bending wave, leaving the structure of the disc at 
the position of the initial pulse relatively unchanged. The 
progression of these two pulses may be seen in the final two 
panels of Fig. (4), where two distinct waves are visible, sep- 
arated by a 'dead-zone' at j/ ~ —0.5. This figure illustrates 
the fact that the code's performance is such that it is able 
to capture the essential features of bending wave propaga- 
tion through Keplerian accretion discs, indicating that the 
vertical structure of the disc models are approximated to 
adequate accuracy. 

In addition to performing bending wave calculations for 
M = 10 discs, a calculation for a TVf = 30 disc was also per- 
formed, with go = 5° in equation (p^). In this case, the value 
of the dimensionless viscosity coefficient a is expected to be 
larger than the disc aspect ratio, since calibration experi- 
ments indicate that a ~ 0.04. These calibration calculations 
measure the value of a acting through the (r, (f)) component 
of the viscosity tensor, rather than that which acts on the 
vertical shear, but should nonetheless provide an estimate 
of the ai acting on this vertical shear. It is expected that a 
warping disturbance in a disc with qi > H/r will evolve dif- 
fusively, rather than in a wave-like manner, and this is what 
is observed. The non linear SPH calculation was fitted using 
the equations from linear theory, and the best fit indicated 
that ai ~ 0.2 for an assumed value of the midplane Mach 
number A4 = 30. There are probably two reasons why this 
somewhat large value is found. 

Firstly the initial pulse is non linear in this case because 
go is about three times larger than H/R ~ in this case, 
whereas it is comparable when M — 10. It thus corresponds 
to go = 15° for the latter Mach number. For such non linear 
warps the indication is that the evolution is more diffusive 
(see below) so resulting in a larger measured value for a^. 
In addition, the vertical structure in the inner parts of the 
disc in this calculation is not as accurately modeled as in 
the AI = 10 calculation, because the smoothing lengths are 
on the order of the disc thickness. This has the effect of 
weakening the pressure forces that arise when the disc mid- 



Tllt Angle v^s Radius 



< 



o 






\ 


0.00 
2.52 




in 








3 47 
4.39 




o 














y — - 

.7 ■ / 






















o 













0.2 0.4 0.6 0.8 1 



Radius 

Figure 7. This figure shows the tilt angle versus radius, as a 
function of time, during a non linear SPH calculation in which 
the initial tilt 30 = 20 degrees. Note the non wave— like evolution 
compared with previous, low amplitude runs. 

planes at adjacent radii are misaligned due to warping, and 
thus reducing the efficacy of both diffusive and wave-like 
communication. This increases the apparent magnitude of 
Qi above that expected from calibration experiments that 
measure the (r, </)) component of viscosity. 

5.2.2 Low Amplitude Tilt in Discs Orbiting in a Softened 
Potential 

A number of calculations were performed for discs which 
were orbiting in softened Keplerian potentials. The discs in 
these calculations initially had an outer radius i? = 1, and 
extended in to i? = 0. The gravitational softening length 
was 6 = 0.2. Here we describe two calculations in which 
M = 10 and the tilt amplitude was go = 7°. One calculation 
employed — 20, 000 particles, whereas the other employed 
N — 102, 000, so that the purpose of this comparison is 
to examine how the results vary as a function of particle 
number. In addition, we are able to observe the effects of 
a small departure from a Keplerian rotation profile induced 
by the gravitational softening. 

The results of the calculation employing N — 102, 000 
particles are plotted in Fig. (5) along with the results from 
a linear calculation. It may be seen from these results the 
substantial effect of changing from a Keplerian potential to 
a softened potential. The linear calculation (represented by 
the dot-dashed line) adopted a Mach number = 10 and 
Qi = 0.04/r^''^. The solid line shows the evolution of g{r) in 
the SPH calculation. The fact that the initial pulse separates 
less cleanly into an inward and outward moving disturbance 
indicates a greater and dominant dispersion resulting from 
the non zero value of — k in this softened potential case. 
The relatively poor fits at radii R > 1.1 at the latest time 
are probably due to the effects of a pressure expansion at 
the outer edge during the SPH calculation, which produced 
a low density region there. As in the previous example of a 
disc orbiting in a Keplerian potential, the drop in the mag- 
nitude g around the position of the initial pulse is not well 
reproduced by the SPH calculation, and this may be again 
due to a localised larger effective viscosity acting in the SPH 
run. Overall, however, the broad features of the solution ob- 
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Figure 8. This figure shows the perturbed radial velocity arising from bending wave propagation through the disc with qq = 20 degrees. 
Note the vertical shear as perturbed radial motion, which is odd in z, is generated in the disc by the warp. The arrows at the bottom of 
each panel indicate the magnitude of the midplane sound speed at each radius. It is apparent that sonic radial motion is generated in 
the disc by the initial large amplitude warp. 



tained using the linearised equations are reproduced by the 
SPH calculation, particularly in the inner parts of the disc, 
indicating that bending wave propagation may be modeled 
using SPH with reasonable accuracy. 

In Fig. (6), we plot the results from the two non linear 
SPH calculations with different numbers of particles. Mov- 
ing through the panels it is apparent that the calculations 
are very similar in terms of the time evolution of the tilt, 
and indicate that accurate results may be obtained with 
particle numbers as low as = 20, 000 for these particular 
calculations. This is because the vertical structure is well 
enough resolved over most of the radial extent of the discs 
for M = 10. 



5.2.3 Evolution of Nonlinear Warps 

In order to examine the effect of increasing the amplitude 
of the initial warping disturbance on the characteristics of 
wave propagation, a number of calculations were performed 
in which the initial warp amplitude was varied. The cal- 
culations were performed with discs orbiting in Keplerian 
potentials, whose midplane Mach number was M — 10, and 
which employed A'^ = 102, 000 p articles (i.e. identical to the 
discs described in section |5.2.1|). The values of go in equa- 
tion @ took the values of 5°, 10°, and 20°. The calcu- 
lation for which go — 5° was described in section (5.2.1). 



Although not described here in detail, the calculation per- 
formed with go = 10° showed behaviour intermediate be- 
tween the go ~ 5° and go — 20° cases, as expected. The evo- 
lution of the tilt angle g{r) as a function of time is displayed 
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Figure 9. This figure shows the full velocity field for the high 
amplitude bending wave calculation. Note the convergence point 
at J/ = —0.5, corresponding to a shock. 



in Fig. (7) for the calculation with go = 20°. Comparing this 
figure with Fig. (1) shows the difference between the go = 5° 
case, in which the warp behaves like a linear perturbation, 
and the case with go = 20°, in which the warp is a non linear 
perturbation. In particular, it may be seen that the evolu- 
tion of a free {i.e. unforced), non linear warping disturbance 
is much less wave-like than for the linear disturbance, in- 
dicating that non linearity in the system leads to enhanced 
dissipation of the bending waves. The evolution of the ve- 
locity field for the calculation in which go = 20° is shown in 
Fig. (8), and should be compared with Fig. (4) which shows 
the same for a linear wave (i.e. go = 5°). It is obvious that 
the strong misalignment of the disc midplane arising from 
the go — 20° warp leads to a flow convergence with radial 
velocities on the order of the sound speed, such that a shock 
forms at r ~ 0.5. This shock does not appear in Fig. (4) for 
which go — 5° , because the degree of warping is insufficient 
to generate large enough radial motions in the disc. The ex- 
istence of this shock, which can be more clearly observed in 
Fig. (9), leads to a strong dissipation of the bending wave, 
such that its temporal evolution appears more diffusive than 
is the case with the go ~ 5° pulse. We note that the existence 
of the shock in this calculation arises because of the form 
of the initial pulse applied to the disc, which automatically 
leads to convergence of the radial flow in the disc. 

In a situation where the warping of the disc is monotonic 
rather than being in the form of a localised pulse, there 
would be no such convergence point, and presumably this 
type of shock would not occur. However, if the degree of 
local warping of the disc is such that radial velocities on the 
order of the sound speed are generated (i.e. when the local 
curvature r\^\ > H/r), then bending waves are unable to 
sustain themselves since the disc material is induced to move 
faster in the radial direction than the wave can propagate. 
It is likely that such a situation will also lead to enhanced 
dissipation in the system, such that the non linearity of the 
warp will be damped. 



In this paper we have performed nonlinear simulations of 
bending waves in accretion discs. Disc Mach numbers rang- 
ing between 10 and 30 have been considered with particle 
numbers ranging between 20000 and 102000. The main pur- 
pose of these calculations is to examine the ability of SPH 
to model warped accretions discs, a problem which is intrin- 
sically three dimensional in nature. 

Linear bending waves are expected to propagate if the a 
viscosity appropriate to dissipation of vertical shear, ai is < 
H/r. In our simulations ai ~ 0.04 for discs with M = 10, so 
bending waves should propagate in this case. We performed 
tests that showed that the linear wave propagation found 
in the simulations could be well fitted by solutions of the 
linearized problem found using a finite difference scheme. 
This is true for strictly Keplerian and slightly non Keplerian 
discs, where the major difference observed between these two 
cases is a more dispersive evolution of the bending waves 
when the central potential is non Keplerian. 

A transition from wave-like to diffusive propagation of 
warps is expected to occur when ai > H/r. Such a transi- 
tion in behaviour was observed in the non linear simulations 
when the disc Mach number was increased to = 30, since 
we estimate that qi > H/r in this case. We therefore re- 
mark that SPH simulations can model warped discs in both 
the wave-like and diffusive regimes. 

A detailed examination of the velocity field produced 
by the simulations showed that the vertical shear motion 
predicted by Papaloizou & Pringle (1983) was reproduced, 
indicating that the vertical structure of the discs was ade- 
quately modeled. 

Calculations of non linear warping disturbances were 
performed, where the condition for non linearity is that the 
local curvature r|^| > H/r. It was observed that non lin- 
ear warps lead to an increased dissipation due to shocks, 
restricting the perturbed horizontal motions to be subsonic. 
This was associated with a transition from wave-like to more 
diffusive-like propagation of the warps, which has the same 
effect as would be produced by increasing the effective value 
of Qi. We speculate that the effects of non linearity will al- 
ways be to inhibit communication between different radial 
locations in the disc which require horizontal velocities that 
are supersonic for the communication to be maintained. We 
expect the effects of non linear damping to adjust the warp 
amplitude until it becomes a linear perturbation. This be- 
haviour is similar to that produced by an increase of the 
effective value of ai. 

The simulations that we have presented in this paper 
indicate that SPH is able to model warped accretion discs 
in the regime where warps are expected to propagate as 
bending waves, and also in the regime where the propaga- 
tion is diffusive. The modeling of warped discs represents a 
stringent test of any hydrodynamical code, since it requires 
the accurate modeling of complex motion in three dimen- 
sions. We find that accurate results may be obtained even 
for rather modest numbers of particles {i.e. N ^ 20000). We 
have previously used SPH simulations to examine the struc- 
ture of warped accretion discs in misaligned binary systems 
(Larwood et al. 1996, Larwood & Papaloizou 1997). We have 
also examined the structure of warped discs orbiting about 
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Kerr black holes (the 'Bardeen-Petterson effect') in a com- 
panion paper to this one (Nelson & Papaloizou 1999). 
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